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Abstract We consider a continuous-time Markov chain (CTMC) whose state 

space is partitioned into aggregates, and each aggregate is assigned a proba- 
bility measure. A sufficient condition for defining a CTMC over the aggregates 
is presented as a variant of weak lumpability, which also characterizes that the 
measure over the original process c;an be recovercsd from that of the aggregated 
one. We show how the applicability of de-aggregation depends on the initial 
distribution. The application section is a major aspect of the article, where we 
illustrate that the stochastic rule-based models for biochemical reaction net- 
works form an important area for usage of the tools developed in the paper. 
For the rule-based models, the construction of the aggregates and computa- 
tion of the distribution over the aggregates are algorithmic. The techniques 
are exemplified in three case studies. 
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Introduction 

The theory of Markov processes has a wide variety of apphcations ranging 
from engineering to biological sciences. In systems biology appropriate Markov 
processes are used in stochastic modeling of different biochemical reaction 
systems, especially where the constituent species are present in low abundance. 
Aggregation or lumping of a Markov chain is instrumental in reducing the size 
of the state space of the chain and in modeling of a partially observable system. 
Typically, the original state space, 5', of the Markov chain {X„} is partitioned 
into a set of equivalence classes, S = {Ai, . . . ,Am}, and a process, {F„}, is 
defined over S. More precisely, let tt be an initial distribution on S for the 
chain {X„}. For a given partition S of S, let the aggregated chain {Yn} be 
defined by 

{Yn = Ara} if and only if {X„ e A™}. 

Observe that {Yn} is not necessarily Markov, nor homogeneous. Conditions 
are imposed on the transition matrix of the Markov chain {X„} to ensure 
that the new process {Yn} is also Markov (see [20], [19], [21], [3], [22] and 
references therein). In this context, strong lumpability refers to the property of 
{X„}, when the aggregated process {Yn} (associated with a given partition) 
is Markov with respect to any initial distribution tt. If P denotes the transi- 
tion matrix of {X„}, then it has been shown that a necessary and sufficient 
condition for {A^n} to be strongly lumpable with respect to the partition S 
is that for every Ak,Ai, J^seAi P{s\s) = J2seA, Pis",s) for any s',s" & Ak- 
Tian and Kannan [55] extended the notion of strong lumpability to continuous 
time Markov chains. A more general situation is when {A^n} is weakly lumpable 
(with respect to a given partition), that is, when {1^} is Markov for a subset 
of initial distributions tt. The notion first appeared in [16] and subsequent 
papers |20 p i9 [ [T7] focussed toward developing an algorithm for characterizing 
the desired set of initial distributions. The characterization is done through 
some kind of recursive equations which sometimes might be hard to read. 

The sufficient condition that we provide in the current paper for {A„} to 
be weakly lumpable with respect to partition S is easy to read and is geared 
toward applications in combinatorial reaction networks. In particular, our con- 
dition enables us to recover information about the original Markov chain from 
the smaller aggregated one (see Theorems [2] [6] [9] [TT]) . This 'invertibility' prop- 
erty is particularly useful for modeling protein networks and is not addressed 
explicitly for weakly lumpable chains in previous literature. A variant of our 
condition can be found in [3] where the author considered backward bisimu- 
lation over a class of weighted automata (finite automata where weights and 
labels are assigned to transitions). For each i, let be a probability mea- 
sure over Ai. The condition that we impose requires that for every i and j, 
X^seA cti{s)P{s, s') / aj{s'), s' e Aj is constant over Aj. The condition can be 
interpreted as follows: Suppose that you are at the state s' € Aj and you look 
back and try to compute the probability that your immediate previous posi- 
tion was somewhere in Ai. The above condition implies that this probability 
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is same no matter where you look back from in Aj. This in particular gener- 
alizes the notion of exact lumpability which corresponds to the case when the 
measures are uniform [3;. Interestingly, if the initial distribution 'tt respects 
a^' in the sense that tt{s) /T:{Ai) ~ ai(s), then the conditional probability 
P{Xt = s \ Yt = Ai) = ai{s), for all t > 0. In fact, we proved that even if 
the initial distribution does not respect the ai, the above result holds asymp- 
totically. These convergence results established in the article are particularly 
useful for modeling purposes and to the best of our knowledge have not been 
discussed before. They imply that the modeler can run the 'smaller', aggre- 
gated process {Yt} and can still extract information about the 'bigger' process 
{Xt} if the need arises. This is further illustrated in the application section. 

The main practical difficulty in aggregating and de-aggregating a Markov 
chain is to construct the appropriate partition, and to find the probability 
measure over the aggregates. Both issues are successfully resolved in the ap- 
plication to the rule-based-models of biochemical reaction networks. 

Traditional modeling of biochemical networks is centered around chemical 
reactions among molecular species and a state of a network is a multi-set of 
molecular species. A species can be, for instance, a protein or its phospho- 
rylated form or a protein complex that consists of several proteins bound to 
each other. Especially, in cellular signal transduction the number of different 
such species can be combinatorially large, due to the rich internal structure 
of proteins and their mutual binding [T3],[53|. For example, one model of the 
early signaling events in epidermal growth factor receptor (EGFR) network, 
with only 8 different proteins gives rise to 2748 different molecular species 
[2]- In such formal description of the cellular process using different 

reactions and species becomes computationally expensive. 

Instead, an efficient way to encode different molecular interactions is to 
use a site-graph based model. A site-graph is a generalization of a graph where 
each node contains different types of sites, and edges can emerge from these 
sites. Molecular species are often suitably represented by site-graphs, where 
nodes are proteins and their sites are the protein binding-domains or modi- 
fiable residues; the edges indicate bonds between proteins. Every species is a 
connected site-graph, and in accordance with the traditional model, a state of 
a network is a multi-set of connected site-graphs. Importantly, more detailed 
description of the species' structure allows to describe interactions locally, 
between parts of molecular species (sometimes refered to as fragments). For 
instance, it can be stated by one rewrite rule, that any species containing a 
protein of type A can have that protein A phosphorylated. In this case, the 
event of phosphorylation of A is independent of the rest of the species' context, 
i.e. A can equally be part of a dimer (complex of two proteins) or of a very 
large protein complex. It is precisely this independence between the molecular 
events we exploit when aggregating states and constructing a suitable aggre- 
gated process. In the present article we present a rigorous construction of a 
Markov chain {Xt} on an appropriate space of site-graphs which essentially 
tells us how the 'reaction soup' looks like at different points of time. It is then 
shown that the usual species-based Markov chain can be constructed as an ag- 
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gregation of {Xt}. But more importantly, there exist other aggregations which 
lead to Markov chains living on much smaller state spaces, and information 
about the species-based model can be extracted at any point of time from 



these smaller Markov chains (see Theorem 12). 

One important feature of the presented application is that it provides an 
effective way of constructing the partition and the accompanying distributions 
over the aggregates (also see [S], [T5]). In particular, the three case studies 
presented at the end exemplify our approach to effectively reduce the state- 
space of the CTMCs in the context of molecular interactions (see Table [l] for 
an overview of the achieved reduction) . 

The work presented in this paper is inspired by the related work of [8], 
where the general algorithm for reducing the stochastic behavior of any Kappa 
[7] model is shown. The proof uses a cumbersome object of a weighted labeled 
transition system, supplied with all the details which are necessary when pro- 
viding a general reduction algorithm. In contrast, in the present article the 
mathematical treatment of the rule-based models has been carried out ef- 
ficiently by using the tools of graph theory. The analysis of Markov chain 
aggregation is done for general Markov chains whose application covers but 
is certainly not limited to the class of rule-based models. In such a set-up, 
the existing reduction framework is extended with a criterion on the rule-set 
for claiming the asymptotic possibility of reconstruction of the species-based 
dynamics. 

The rest of the work is organized as follows. In Section [T] we describe con- 
ditions on the transition matrix and initial distribution of the Markov chain 
{Xn} which will ensure that the aggregated chain {Yn} is also Markov. The 
conditions described are tailor-made for our applications to biochemical reac- 
tion networks. We also prove convergence properties of the transition proba- 
bilities of the aggregated chain when the initial distribution does not satisfy 
the required conditions. The case of continuous time chains has been treated 
in Section [2] Section [3] first discusses the traditional Markov chain modeling 
of biochemical reaction systems using reactions and species. Next, the math- 
ematical definition of site-graphs is introduced and the formal description of 
site-graph based modeling of protein-protein interaction is given. Section [4] is 
devoted to applications. We describe the criteria for testing the aggregation 
conditions on the CTMCs which underly rule-based models. Illustrative case 
studies are given at the end. 



1 Discrete time case 



Let {^n} be a Markov chain taking values in a finite set S with transition 
matrix P and initial probability distribution tt. Let S = {Ai, . . . ^ Am} be 
a finite partition of S. Moreover, let {ai}i=i...._m be a family of probability 
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measures on S, such that ai{s) — ioi s ^ Ai. Define S : S x S ^ M>o by 

6{Ai, s) = '- r—^ , where s d An. 

aj{s) 

Assume that the following condition holds. 

(Condi) For any A;, Aj e S and s, s' e Aj, S{Ai, s) = S{Ai, s'). 

Fix s G Aj and let P{Ai,Aj) := S{Ai,s). Notice that P is unambiguously 
defined under [(Condi) 



Theorem 1 P is a probability transition matrix. 



Proof Notice that by (Condi) 



A,{s)P{A,,A,)^ ^ a,{s')P{s',s). 



Summing over s G Aj, we have 

p{A,A,)^ E 

s'eAi seAj 

It follows that 

rn m 
j = l j = l s'£Ai seAj 

- E a.(^')E^(^''^)= E "»(^') 

s'€Ai seS s'&Ai 

= 1. 

Definition 1 For any probability distribution tt on S, define the probability 
distributions 7r|^. on Ai and tt on S" by 



Definition 2 We say that a probability distribution tt respects {ai : i 

1, . . . , m} if tt\j\. (s) = Q!i(s) for s £ Ai, i — 1, . . . ,m. 



1.1 Aggregation and de-aggregation 

Throughout, we will assume that {Yn} is a Markov chain taking values in S 
with transition matrix P and initial distribution tt. 

Theorem 2 Assume that tt respects {ai : i = 1, . . . ,m}. Then for all n — 
0,1,.-. 
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(i) (lumpability) P(y„ = A,) = P(X„ e A,); 

(ii) (invertibility) P(X„ = s) = P(r„ = A,)a,{s). 



We need the foUowing two lemmas to prove Theorem [2] 



Lemma 1 Assume that for all i = l,...,m, P(X„_i = s|X„_i <E Ai 
^P""iU.(s) = a,is). Then P(X„ e A,|X„_i e A,) = P(A„ A,). 



Proof Notice that 



P(x„eA,|x„_ieA,) = 



J2s'eA, J2seA, ^i^n - S,X„-1 - s') 

P(X„_i e A,) 

P(x„_i e A,) 

E.^gA. EsgA, P(^»-i e A,)P(X„_i ^ g'|X„_i e A,)P{s\s) 

P(X„_i G AO 

P(x„_i e AO 



by the hypothesis 



= P(Ai, Aj)aj(s), by the definition of P 

seAj 

= P(A„A,) J2 aj{s) = PiA,A,). 
seAj 



Lemma 2 Assume that 7r|^. (s) — ai{s) for s e A,;,i = l,...,m. Then 



Proof The first equality is of course by the definition. For the second, we 
use induction. The case n = is given. Suppose that the statement holds for 
k — n—1. First observe that if s ^ A^, then both sides equal 0. So assume that 
s e Ai. Then, by Lemma[l] we have that P(X„ e Aj\Xn-i e AO = P{Ai,Aj). 
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Next note that 



P(X„ = 3) 

P(^„ e A,) 

Er=l E.^gA, P(^n-1 € ^,)P(X»-1 = £ A,)P{3\S) 

E;=i P(^n-i e 



ai(s) 



ai(s) 



Er=i P(^n-i e ^j)P(^n e A,|X„_i e A,) 



Q!i(s). 



P(^„ e Ai) 
We next proceed to prove Theorem |2] 

Proof (Theorem [2| We use induction. Notice that both the statements hold 
for n = 0. Assume that (i) and |(ii)| hold for n — 1. Then P(X,i_i — s\Xn-\ G 
A^) = ai{s), and hence by Lemma [l]P(X„ e Aj\Xn-\ e Ai) = P(Ai,Aj). 
Therefore, 



p(r„ = A,) = ^ P(r„_i A,)p{A,,Ai) 

m 

= ^ P(X„_i e yl,)P(X„ e A,|X„_i e A,) 

= P(X„ e A,). 
This proves (i) Next, notice that Lemma |2] implies 
P(X„ = s)=a,(s)P(X„e AO 



a,;(s)P(r„-A,), by (i) 



This proves (ii) 



Remark 1 Notice that we have proved that under the assumption 7r|^;(s) 
a^{s), P(X„ G Aj|X„_i G A,) =P(A,,Aj), for n 1, 2, . . .. 
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1.2 Convergence 

In the previous section, we proved that if {^n} is a discrete time Markov chain 
on S with initial distribution tt respecting {ai : i — I, . . . , m}, then the aggre- 
gate process {1^} is an aggregated Markov chain satisfying lumpabiUty and 
invertibihty property. We now investigate the case when the initial distribu- 
tion of {Xn} doesn't respect {a^ :{ = !,..., m}. We start with the following 
theorem. 

Theorem 3 P"(Ai,Ai) = \ ' for any s S A,. 

Proof We use induction. Notice that for n = 1, the assertion is true by the 
definition of P. Assume that the statement holds for some n. Then, 



P"+\A,A,) J2 PiAr,Ak)P'\Ak, A,) 



EE 



a,{s')Pis',so) 



afe(so) 



E 



akis")P^is",s) 
aj{s) 



E E E 
E E (e 



a,(s')-P(s',so)\ fak{s")P"{s",s) 



aj{s) 



a,is')P{s',s")\ (a k{s")P^{s\s) 



^E E fE"^(-w,^'V"(^",^)) 



for any sq G Ak 



by (Condi) 



aAs, 

■' ^ ' s'^Ai k s"eAk 

aj{s) 

We say that s — s' , if for some n > 0, P"(s, s') > 0. Recall that the Markov 
chain {X„} is irreducible if s s' for any s,s' G S. One corollary of Theorem 
[3] is that if for s G Ai, s' & Aj , s ^ s' then Ai -> Aj for the Markov chain Y. 
In fact, we have the following result. 

Theorem 4 Let {Xn} be a discrete time Markov chain on S with transi- 
tion probability matrix P and {Yn} a Markov chain taking values in S with 
transition matrix P. Then 

(i) If the process {Xn} is irreducible, then so is {Yn}. 

(ii) If s S is recurrent for the process {Xn}, then so is Ai for the process 

{Yn}- 

(iii) li s E Ai has period 1, then the period of Ai is also 1. 
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Proof (i) is immediate from the discussion. Notice that if s E Ai,s' e Aj, 
then by Theorem [sj^ P" (A, , ) > ^j^P"(s,s'). Therefore, it follows that if 
s € Ai — Aj, then P''^{Ai,Ai) > P"{s, s). Now s € Ai is recurrent if and only 
if J2n ^) — Observe that 



It follows that for the process {Yn}, Ai is recurrent. This proves (ii) Next 
observe that 

{n : P"(s,s) > 0} C {n : P'\A,A,) > 0}, 



and (iii) follows immediately. 

For the following results we will assume that there exists a probability 
distribution vr on 5 which respects {cti : i = 1, . . . , m}. 

Theorem 5 Let {Xn} be a discrete time Markov chain on S with transition 
probability matrix P and unique stationary distribution /i. Then /i respects 
{a,; : i = I, . . . , m}. 

Proof Let tt be a probability distribution on S which respects {ai : i — 
1, . . . , m}. Now since /i is unique, we have for any set A, ^ X]I-=i 7rP'^(A) 
u{A) (see [13]). By the choice of tt, 7r(s) = ai{s)TT{Ai) for s G Ai. By Theorem 
|2j t:P''{s) = a^{s)TTP''{A,), s € A,. Therefore, it follows that ^ X^Li 7rP'=(s) = 
ai{s)^ T^P'^iAi), s E Ai. Taking limit as n cx), it implies that /ij^i; (s) = 

ai(s), s € Ai,i ^ 1, . . . ,m. 

For any set A C S, let P(")(s,A) := - P'' {s , A) , and P(")(^n S 

Theorem 6 Let {Xn} be an irreducible Markov chain taking values in S 
with transition matrix P. Let fj, be the stationary distribution of P. Let {Yn} 
be a Markov chain on S with transition matrix P. Then /i is the stationary 
distribution for P. Also for all n = 0, 1, . . ., 

(i) p(")(r„ = A,)-p("n^ne^,)^0; 

(ii) P(")(X„ = S)/P(")(r„ = A,) ^ a,(s). 

Proof We first show that /i is a stationary distribution of P. Towards this end, 
first observe that by Theorem [5] /i respects {ai : i = 1, . . . , m}. Take /i as the 
initial distribution of {Xn}. Then by[(i)]of Theorem [2] 

m(A) = fi{A,) = /iP"(A) = M^"(AO. 

It follows that p, is stationary for P. Since P is irreducible by Theorem |4j jl 
is unique. Now let tt be any initial distribution for P. S ince i s th e unique 
stationary distribution for P, 7rP("^(Ai) -> /i(Aj). Hence (i) and (ii) follow. 



The above result can be improved if we assume in addition that the Markov 
chain {^n} is aperiodic. 
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Theorem 7 Let {X„} be an irreducible, aperiodic Markov chain taking values 
in S with transition matrix P. Let be the stationary distribution of P. Let 
{Y-n} be a Markov chain on S with transition matrix P. Then 

(i) P(y„ - A,) - P(X„ e A,) ^ 0; 

(ii) P(X„ = s)/P(r„=yl,)^a,(s). 

Proof By Theorem]!] the Markov chain {Yn} is also aperiodic and irreducible. 
Moreover by the previous theorem, // is the unique stationary distribution for 
{y„}. The result follows by noting that for any aperiodic, irreducible Markov 
chain {Zn} with a stationary distribution rj, P{Zn G A) — > vi^)- 

2 Continuous time case 

We now consider a continuous time Markov chain, {^t}fg[o,oo)j taking values 
in a countable set S. Let Q be the generator matrix for {Xt}. As before, let 
S = {Al, . . . , A^} be a finite partition of S and {ai} be a family of probability 
measures on S with ai{s) = O,for s ^ Ai. Define A : S x S ^ K>o by 

s) = ^'^^^ \ , ^ \ where s G A,. (1) 
Assume the following condition holds. 

(Cond2) For any Ai,Aj e S and s, s' e Aj, A{Ai, s) = A{Ai, s'). 

Fix s e and let Q{Ai,Aj) := Z\(Ai,s). Notice that Q is unambiguously 
defined under [(Condi) 

Theorem 8 Q is a generator matrix. 

Proof We only need to prove that X^Jli Q{^i:^j) = 0. The proof proceeds 
almost exactly in the same way as that of Theorem [T] 

For any generator matrix Q = (qij), define 

3 

2.1 Aggregation and de-aggregation 

We next prove the analogue of Theorem [2] 

Theorem 9 Let {Xt] be a continuous time Markov chain taking values in a 
countable set S with generator matrix Q and initial probability distribution tt. 
Let {Yt \ be a continuous time Markov chain taking values in S with generator 
matrix Q and initial distribution tt. Assume that tt respects {oii : i = 1, . . . , m}. 
Also assume that there exists an r > such that supj qi < r. Then for alH > 
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(i) (lumpability) P{Yt = A,) = P{Xt G A,); 

(ii) (invertibility) P{Xt = s) = P{Yt = Ai)a^{s). 

We prove the above theorem by constructing a uniformized discrete time 
Markov chain out of {Xt}. For any matrix A = ((aij))^ jgs, we use the norm 
1 1 ^11 = supj^^ I fly I . Note by the assumptions in Theorem |9j \\Q\\ < r < oo. 
If P denotes the transition probabihty matrix of {Xf}, then P satisfies the 
Kolmogorov forward equation 

P'{t) = P{t)Q,t> 0. 

Since \\Q\\ < oo, the solution to the above equation is given by 

oo 

p{t) = = Y^m'/ki 

k=0 

Define the transition matrix M by M = I + Q/r. Writing Q — r{M — I) we 
have 

fc=0 

Let {Zn} be a Markov chain on S with transition probabihty matrix M. Let 
^ be a Poisson process with intensity r independent of {Zn}- Then ([2| imphes 

that {Xt} = {Z{^{t))}. We wih need to consider the aggregate Markov chain 
{Zn} on S with the transition matrix defined by 



M{A,A,) = ) ^ ' \ se A,. (3) 

0:j(S) 



Lemma 3 M is well-defined. 

Proof We need to show that M{Ai,Aj) does not depend on the choice of 
s e Aj. Let p{Ai,s) denote the right side of ([3| for s g Aj. We will use 



(Cond2) First assume that i ^ j- Then, 

PiA,s) = l ^^'eAM^Qi^'^s) ^ lQ^A.,A,),se A,. (4) 



by the definition of the Q matrix. If i = j, then 

P[A^,S) = ^- 

ai(s) 

_ E.'sA. a^{s')Q{s\ s)/r - a^{s)Q{s, s)/r + a,{s){l + Q(s, s)/r) 
1 y^»'c4 ai{s')Q{s's) 1 . 

It follows p{Ai^ s) is independent of the choice of s G A^, j = 1, . . . , m. 
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We now prove the foUowing commutativity relation. 

Theorem 10 Let Q be a generator matrix with sup^ qi < r for some r, and 
let 

(i) {Xt} be a continuous time Markov chain taking values in a countable 
set S with generator matrix Q and initial probability distribution tt. 

(ii) {Yt} be a continuous time Markov chain taking values in S with gener- 
ator matrix Q and initial distribution tt. 

(iii) {Zn} be the uniformized discrete time Markov chain (corresponding to 
{Xt}) on S with transition matrix M = I + Q/r and initial distribution 

TT. 

(iv) {Yn} be the uniformized discrete time chain (corresponding to {5^t}) on 
S with transition matrix M = I + Q/r and initial distribution tt. 

(v) {Zn} be the discrete time Markov chain on S with transition matrix M 
and initial distribution tt. 

Then {ij = {y„}. 

Proof We only need to show that M = M . But this readily follows from Q 
and 

Proof (Theorem |9| Note that ^ implies. 



k\ 

i:p,z....)q?£ 

fe>0 



e-''\rt)'' 



1:p(Zi6A,) 

k>a 



seAi k>a 



PiXt = s) = PiXt e A, 



seAi 



Here, the second equality is by Theorem 10 while the third is by (i) of Theorem 
[2) This proves [(i)] and | (ii) | follows similarly. 



2.2 Convergence 

Let /i be a stationary distribution of the continuous time Markov chain {Xt}, 
that is /i satisfies /zQ = 0. Then we have the corresponding analogue of The- 
orem |6l 
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Theorem 11 Let {Xt} be an irreducible Markov chain taking values in S 
with generator matrix Q. Assume that sup, Qi < r, for some r > 0. Let fi 
be the stationary distribution of Q. Let {Yf} be a Markov chain on S with 
generator matrix Q. Then fl is the stationary distribution for Q. Moreover, 

(i) P{Yt = A,) - PiXt e A,) ^ 0; 

(ii) P{Xt^s)/PiYt=A,)~^a,is). 

Proof We first consider the uniformized chain {Zn} corresponding to {Xt} 
with transition matrix AI = I + Q/r. Note that fi is the stationary distribu- 
tion for {M}. It follows by Theorem [6j that p, is the stationary distribution 
for {Zn}, hence for {F„}- It follows that flQ = 0. Next supj Qi < oo guarantees 
that the chain does not explode. The result follows by noting that for any irre- 
ducible, non-exploding continuous time Markov chain {Zt} with a stationary 
distribution 77, P{Zt G A) — >■ t]{A) as t — > 00. 

3 Formalism 

The standard model of biochemical networks is typically based on counting 
chemical species (complexes). However, for our purpose it is useful to consider 
a site-graph based description of the model. We start by briefly outlining the 
Markov chain formulation of a species-based model of a biochemical reaction 
system, and then move on to the concept of site-graph. 

3.1 Modeling biochemical networks by a CTMC 

A biochemical reaction system involves multiple chemical reactions and several 
species. In general, chemical reactions in single cells occur far from thermo- 
dynamic equilibrium and the number of molecules of chemical species is often 
low |15| . |llj . Recent advances in real-time single cell imaging, micro-fluidic 
techniques and synthetic biology have testified to the random nature of gene 
expression and protein abundance in single cells [2S] , [5] . Thus a stochastic de- 
scription of chemical reactions is often mandatory to analyze the behavior of 
the system. The dynamics of the system is typically modeled by a continuous- 
time Markov chain (CTMC) with the state being the number of molecules 
of each species. [Tj is a good reference for a review of the tools of Markov 
processes used in the reaction network systems. 

Consider a biochemical reaction system consisting of n species and v reac- 
tions, and let X{t) denote the state of the system at time t in Z" . If the fc-th re- 
action occurs at time t, then the system is updated as X{t) = X{t—) + i''^ — v^, 
where X{t—) denotes the state of the system just before time t, and i'^ ^ 
Z" represent the vector of number of molecules consumed and created in one 
occurrence of reaction /c, respectively. For convenience, let vi^ — v'^ ^ ■ The 
evolution of the process X is modeled by 



P[X{t + At)^x + Vk\X{t) = x] = ak{x)At + o{At). 
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The quantity ak is usually called the propensity of the reaction k in the chem- 
ical literature, and its expression is often calculated by using the law of mass 
action [21], [ID]- The generator matrix or the Q-matrix of the CTMC X is 
given by qx,x+Vf, = o,k{x). The CTMC X will have an invariant measure tt if 
ttQ = 0. 



3.2 Site-graphs 

The notion of a site-graph is a generalization of that of a standard graph. A 
site-graph consists of nodes and edges; Each node is assigned a set of sites, and 
the edges are established between two sites of (different) nodes. The nodes of 
a site-graph can be interpreted as protein names, and sites of a node stand for 
protein binding domains. Let S denote the set of all the sites in a site-graph, 
and let V{S) denote the the class of all subsets of S. 

Definition 3 A site-graph G = (F, i7, E) is defined by a set of nodes V", an 
interface function U : V ^ and a set of edges E C {{{v, s), {v' , s')}\v, v' S 

V,v:^v',se IJ{v),s' e U{v')}. 

The function U in the above definition tracks the sites corresponding to a 
particular node of a site-graph. 

Definition 4 Given a site-graph G = (V, Z", a sequence of edges (ei, ... e^) £ 
E'', Ci = {{vi, Si), {v'^, s'i)}, such that v'^ = w^+i and s'^ ^ s^+i for i = 1, . . . fc— 1, 
is called a path between nodes v\ and f fc . If there exists a path between every 
two nodes v,v' €V, a site-graph G = {V,S,E) is connected. 

Definition 5 Let G = {V, S, E) be a site-graph. A site graph G' is a sub- 
site-graph of G, written G' C G, if V C V, for all v G V , E'{v) C i:{v), 
and E' C E. 



3.3 Site-graph-rewrite rules 

Definition 6 Let G = (V, E) be a site-graph. We introduce two elementary 
site-graph transformations: adding/deleting an edge. 

*^ae(G, 6) — (Vnew T ^ : -^new ) • ^new — ^new — EU {c}, 
^de(.G, e) — {^new 7 -Enevo)- ^new — F; -^new — \ 

The interface function S is unaltered under any of the above transformations. 
Let G' = {V',S,E') be a site-graph derived from G = {V,S,E) by a finite 
number of applications of S^ni <^ae, Sde- Let c G K>o be a non-negative real 
number denoting the rate of the transformation. The triple (G, G',c), also 
denoted by G A G', is called a site-graph-rewrite rule. 
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3.4 Rule-based model 

Suppose that TZ = {Ri, . . . , i?„} is a collection of site-graph rewrite rules such 
that for i — 1, . . . ,n, Ri = (Gi, G^, Ci) and Gi = {Vi, Si, Ei). From now on, for 
a given set of rules TZ, we use the terminology 

— the set of node types for V U^Vi, 

— the set of edge types for E :~ UiEi, 

— the interface function for S : V ^^{S), such that for v £ V , S{v) := 

For each node v G V, we will consider Uy copies or instances of the node v, 
denoted hy ,v'^ , . . . , w"^' . Note that, in the Kappa rule-based models, the set 
of node types and edge types are predefined in the signature of the model; 
Here, it is deduced from the set of rules (a more detailed discussion to the 



relation with Kappa is given in Section 5.51 



Definition 7 A reaction mixture is a site-graph Q — {V,S,£) where 

- V^{v^veV,j = l,...,ny}; 

- £ C {{(«!, si), (w^,S2)}|{(wi,si), {v2,S2)} e E,i = l,...ny,,j = l,...,ny^} 

Definition 8 A rule-based model is a collection of rules TZ, accompanied 
with the initial reaction mixture Qq. 

Remark 2 By definition, the site-graphs Gi and occurring in some rule 
(Gi, G'i,Ci), are such that a node v £ V, edge e £ E, but also a site s G S{v) 
may be omitted: for some rule Ri, we may have a node v € Vi, such that there 
exists a site s G S{v) \ Si{v). The possibility of omitting a site s € S{v) from 
the interface of node v means that the value of site s does not make an influence 
on the applicability of this rule. This is the crucial aspect of reductions of site- 
graph-rewrite models, because it will help to detect and prove symmetries in 
the underlying CTMC before considering its full generator matrix. 

Definition 9 A rule {Gi, G'i,Ci) is reversible, if there exists a rule {Gj, G'^Cj), 
such that Gi = G'j and G^ — Gj- A rule-based model is reversible, if all its 
rules are reversible. 

Let G be the set of all reaction mixtures which can be reached by finite 
number of applications of rules from 7?, to a reaction mixture Qo- We will now 
describe a Markov chain taking values in Q. The following notion of renaming 
a site-graph will be used for the formal description. 

Definition 10 Let G — {V,S,E) be a site-graph, V' a set such that \V'\ > 
\V\ (I • I denotes the set cardinality), and r/ : V V a,n injective func- 
tion. Then the 77-induced node-renamed site-graph, G^', is given by G^ = 
{viV),U\E^), where S^{r^{v)) = S{v) and E'^ = {{{tjM, s,), (r^iv^), S2)} \ 

vi,v2 e V}. 
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Fig. 1 Case study 1: Simple scaffold, a) The model consists of two reversible rules: a scaffold 
B has two binding sites, a and c, which serve for binding nodes A and C, respectively, b) 
The application of rule _Ri to the reaction mixture Q via node renaming funcion rj results 
in a reaction mixture Q' , which is equivalent to Q except in the sub-site-graph captured by 
node renaming r). 



3.5 The CTMC of a rule-based model 

Consider a reaction mixture G G, a rule Ri — {Gi, G[,Ci) G TZ. Suppose that 
77 : — V is a node renaming function such that C Q. This implies that 
the rule Ri can be applied to a part of the reaction mixture Q. Let Q[-^ ^ be the 
unique reaction mixture obtained after the application of the rule Ri. (For a 
more formal definition of Q'^ , see [6].) Note that G'i C Q' . Define the transition 
rate Q by Q{G, Q'r^^i) = Ci. More precisely, 

{Ci if Q' = Q'^ i for some 77, i 

if - g'^^, for any 77 and i (6) 

~Yjg'^gQiQ^Q') otherwise. 

Let {Xt} be a CTMC with state-space G and generator matrix Q. 
3.5.1 Case study 1: Simple scaffold. 

Consider a site-graph-rewrite model TZ = i?2, -^4)7 depicted in Fig- 
ure We have that V = U^^^^i = {A,B,C}, S{A) = {6}, S{B) = {a,c}, 
SiC) = {6}, and E = {{{A,h),{B,a)},{{G,b),{B,c)}}. In Figure [i]d, we 
show the application of rule Ri to the reaction mixture Q = (V, S, £) , such 
that V = {A\B\B^,B^,C^}, and £ = {{{B^ c) , {C\ b))} , E{A'^) = {&}, 
S{B^) = S{B^) = S{B^) = {a,c}, E{C^) = {b}. 



4 Application 

This section is devoted to establishing applicability of the results from Sec- 
tion [2] to rule-based models. Each of the properties - lumpability, invertability 
and convergence are illustrated on three case studies. For each case study, we 
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(Cl+C2)l2 



Fig. 2 Illustration for testing 



(Cond3) 



and its relation to 



(Cond2) 



Let 



{g,g'}. For g,g' e A2, the permutation a(gi_) = g{, a{g[) = g 



''■{62) = 62 S'nd f (62) ~ ^2 proves that the predecessors of g and those of g' inside class 



Ai are in bijection, that is, (CondS) 
<5(A2,5')> which is because Q(Gi,G) 



holds. 



(Cond2) | follows, since Q(Ai,A2) = <5(Ai,g) = 

Q(<T(gi),e') + Q(o-(62),e') = ci + ca. 



and the rate in the aggregated chain is Q(Ai, Aj) = |^^(ci + C2) = §(ci + C2) 



first define a trivial uniform aggregation of Xt, denoted by Yt, which corre- 
sponds to the usual population-based description with mass-action kinetics. 
We then show that there exists another uniform aggregation of Xt, denoted 
by Zt, with much smaller state space. Finally, since the standard biological 
analysis are referring to the population-based Markov chain we outline below a 
method of retrieving the conditional distribution of Yt given Zt . The summary 
of all considered reductions is given in Table [T] 

The following observation establishes an algorithmic criterion for checking 
(Cond2) and is obvious from ([7|. An illustration is given in Figure [2] 



Lemma 4 Let G = {Ai, . . . , A„} be a partitioning of G induced by an equiv- 
alence relation G x G. Let ai be the uniform probability measure on A^, 
that is, for any Q e A^, ai{Q) = . Note that in this case ([ij reduces to, 

Z\(A,;,g) = |^ ^ Q(C?i,g), g^K,. (7) 



Then, the following condition implies (Cond2) 

(CondS) For all Ai,Aj e G, for all G,G' G Aj, there exists a permutation of 
states in A^, cr : A^ — > A^, such that Q{Qi,Q) = Q{(j{Qi),Q'). 



Definition 11 If the equivalence relation G x G satisfies (CondS) and for 
each i = 1, . . . , m, ai is a uniform probability measure on A^, then the corre- 
sponding Markov chain {Yt} (with generator matrix Q{Ai, Aj) = A{Ai, G),G ^ 
Aj) is a uniform aggregation oi {Xt}. 

Let and ^2 be two equivalence relations of G, such that Gi — {Ai,A2,...} 
and G2 = {Bi,B2, . . .} are the corresponding sets of equivalence classes. Sup- 
pose that ~i and ^2 induce uniform aggregations on {X^}, denoted respec- 
tively by {Yt} and {Zt}. The property of invertibility allows to evaluate the 
conditional distributions of Xt given Yt, and of Xt given Zt. However, as men- 
tioned oftentimes it is of interest to the modeler to retrieve the conditional 
distribution of Yt given Zt. This is possible by the following result. 
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Theorem 12 If is coarser than ^-^2 (that is, ~iC^2), then G2 can be 
obtained by partitioning Gi as follows. 

A, ^ Aj iff there exist G £ A^,g' e Aj, such that Q ^2 G' ■ 

Equivalently, 

Al ^ Aj iff there exists Bfe such that U C . (8) 

Assume that {Yt} and {Zt} with generator matrices Qi and Q2 are two 
uniform aggregations of the Markov chain {Xt} induced by {^i,{ai}) and 
(^2, {A}), where and f3j are uniform over A^ and Mj respectively. Define 

, flM, if A, C B, 
a(A,):= (9) 
, otherwise. 



Then {a'j} satisfies (Cond2) and hence {Yt} is an aggregation of the Markov 
chain {Zt}. 

Proof It is trivial to check that ~ defined by (|8| is a well-defined equivalence 
relation. 

Assume now that Bj,Bj' g G2 and A^/ C Mji. We have to show that 
lj,Air) is constant for all A,/ C Mj/. Toward this end notice that 

, , E^:A.cB,«j(AOQi(A,,A,0 E,:A.cbJA,|/|B,|Qi(A„A,0 



a, (A,,) |A,'|/|B,'| 
EA..cBjA.|/|B,|Eg.gA.Q(g',g)|A.n/|A 

I Ad/1% I 
EA.cB,Es'eA.Q(^^',e)|A. 



for some Q € A^ 



|A,,|/|B,v 
J2 Q{g\5)\M,,\/\M,\ 



6'eBj 
=(52(Bj, Bj'). 

Here the third and the last equalities are because by the assumption {Yt} 
and {Zt} are uniform aggregations of {Xt}. 



4.1 Case study 1: Simple scaffold (continued) 

The simple scaffold example serves as an illustrative case study which demon- 
strates all the introduced concepts in detail. 

Species. A molecular species is a class of connected reaction mixtures which 
are isomorphic up to renaming of the nodes of same type. We here omit 
a formal definition of species, since it is not necessary for conveying the 
arguments. In the scaffold example, all species can be categorized into six 
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Si M5i) MGi) 

Fig. 3 Case study 1: Simple scaffold. A reaction mixture 5i(left) and its graphical repre- 
sentation in aggregation (f>i (center) or 1^2 (right). 



types: (A)- a free node of type A, (B)- a free node of type B, (C)- a free 
node of type C, {AB)~ a node of type B that is bound to node of a type 
A, and is not bound to a node of type C, (BC)- a node of type B that 
is bound to a node of type C, and is not bound to a node of type A, and 
(ABC)- a node of type B that is bound to a node of type A, and is also 
bound to a node of type C. All reaction mixtures G, which count the same 
number of each of the species correspond to the same population-based 
state. The population-based encoding of the state space is captured by the 
function (j>i : G — >■ N^, such that ipiiG) — {'^ab,i^bc-,'>t^abc)i if G has 
mAB sub-site-graphs of type {AB), uibc sub-site-graphs of type (BC), and 
iTLABC sub-site-graphs of type (ABC). Note that, given the value (j)i{G), 
the number of sub-site-graphs of type (A), (B) and (C) in G is also known, 
since the total number of nodes of each type is conserved. Two reaction 
mixtures G and G' are aggregated by relation G x G if they have the 
same value of function 0i : 

G-i G' iff 'Pi{G)^'Pi{G'). 
For example, in Figure [t] (piiGi) 7^ (j)i{G2)- 

The aggregated CTMC, {Yt}, takes values in N"^, and it is exactly the 
standard population-based model description with mass-action kinetics. 
Fragments. The sites a and c of nodes of type B are updated without testing 
each-other. As formally shown later in Lemma[5j any two states which have 
the same number of free sites c and free sites a are not distinguishable by 
the system's dynamics. As a consequence, the following lumping is also 
applicable: let (/)2 : G — )■ be such that 4'2{G) — {mAB*,rn^Bc), if G G 
has ruAB* nodes B bound to A and m^BC nodes B that bound to C. The 
two states G and G' are aggregated by relation G x G if they have the 
same value of function 02 ■ 

G^2 G' iff MG)^<p2{G'). 

For example, in Fig. [jj (j)2{Gi) = 02(^2)- The aggregation of {XJ by 02 
results in a CTMC {Zt], which takes values in N^, and it therefore provides 
a better reduction than the standard population-based model description. 
A way to visualize the states of CTMC's {Xt], {Yt} and {Zt] is shown in 
Fig. [31 
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c 

Fig. 4 Case study 2; Two-sided polymerization. 



Lemma 5 Both relations ~i and ~2 induce uniform aggregations of {^t}. 
Moreover, ^iC^2, that is, ^2 is coarser than ^1. 

Proof Consider lumping by ^2- Let Gi,G2 be two reaction mixtures such that 
Gi ~2 G2, and let 02(Si) = 02(^2) = {rnAB*,rn^,Bc)- If Gi,G2 € Bj, by 
Theorem [12] it is enough to show that for any Bj G G2, and any € B^, 
there is a permutation cr : Bj — > Bj, such that Q{Gi Gi) = Q{(^{G),G2)- Choose 
some Bj e G2 and G G Bj. Then, (j)2iG) e {("i^s* - 1, to*bc), (was* + 
l,"i*Bc), (was*, "i*BC + l), (rriAB*, »7i*sc-l)}- We analyze the case 02 (^) = 
{niAB* — l,w*Bc); the other three cases are analogous. 

Let Gi = {V,S,£i) andfJ2 = {V,S,£2)- Since 02(^1) = 4>2{G2), there exists 
a bijective renaming function ry : V — > V, such that G2 — Gi, that is, G2 is 
induced node-renamed site-graph Gi- It is easy to inspect that Q{GtGi) = ci 
if and only if Q{G^\G2) = ci. So - the bijection over the reaction mixtures 
aggregated to Bj is the one induced by renaming 77. 

For showing that '^2 is coarser than ~i, it is enough to observe that the map 
: -J' defined by (t>{mAB,'rnBc, niABc) = {niAB+mABCi ruBc+mABc) 
is such that (j)2 — 4> ° 4'i- 

Consequently to Lemmajsj Theorem 12 applies. Then, the process {Zt\ is 
also lumpable with respect to {Yt}, and Theorem [9] applies. Imagine that it is 
possible to experimentally synthesize only the complexes of type {AB) and of 
(JSC), but not a complex of type (A), (B), (C) or [ABC). Then, the initial 
distribution does not respect aj, as soon as ua > 1, riB > 2, nc > 1- However, 
since each reversible rule-based model trivially has an irreducible CTMC, the 
Theorem [TT] holds. 

A concrete example is demonstrated in Figure [7] The details for the cal- 
culation for Table [1] de-aggregation, as well the discussion for ua = nc = 1, 
ub = '2- can be found in the Appendix. 



4.2 Case study 2: Two-sided polymerization 

The two-sided polymerization case study illustrates the drastic advantage 
of using the fragment-based CTMC, because it shows to have exponentially 
smaller state space than the species-based CTMC. 

Consider a site-graph-rewrite model TZ depicted in Fig. [SJd: proteins A and 
B can polymerize by forming bonds of two kinds: between site h of protein A 
and site a of protein B, or between site r of protein A and site I of protein 
B. Assume that there are ua nodes of type A and ub nodes of type B. Let 
G be the set of all reaction mixtures. All connected site-graphs occurring in a 
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EGFR 



Fig. 5 a) Summary of interactions between nodes in case study 3. The dotted lines represent 
phosphorylation, and solid lines denote standard bindings. The self-loop at the site d of node 
EGFR means that it can bind to another node EGFR, i.e. receptor dimerization. b) An 
example of a Kappa rule, and a corresponding site-graph rewrite rule. 

reaction mixture can be categorized into two types: chains and rings. Chains 
are the connected site-graphs having two free sites, and rings are those having 
no free sites. We say that a chain or a ring is of length i if it has i bonds in 
total. Chains can be classified into four different kinds, depending on which 
sites are free. 

Species. Let 0i : G ^ N'""™ be such that 
4>i{Q) = (^11,..., a; 
if ^ G G has 

— xii chains of type {A..B)i, that is, of length 2i — 1, with free sites b and 
a, 

— X2i chains of type {B..A)i, that is, of length 2i — 1, with free sites I and 

— x^i chains of type {A..A)i, that is, of length 2i, with free sites b and a, 

— Xi^i chains of type {B..B)i, that is, of length 2i, with free sites I and r, 

— x^i rings of type {.A..B.)i, that is, of length 2i. 

The two states Q and Q are aggregated by the equivalence relation ^iC 
SxSii MG) = MG)- 
Fragments. Let ^2 : G — >■ be such that (j)2{G) — {mri,mba), if ^/ G <G 
has rriri bonds between sites r and I, and m^a bonds between sites b and 
a. The two states G and G' are aggregated by the equivalence relation 

~2eGxGif MG) = MG')- 

Alternatively, since the rates of forming and releasing bonds do not depend 
on the type of the bond, let 1^3 : G — >■ N be such that 03(t/) = m, if C/ G G 
has in total m bonds. The two states G and G' be aggregated by equivalence 
relation -36 G x G if MG) = MG')- 

A concrete example is demonstrated in Figure [8] The details for the calcu- 
lation for Table [U and on de-aggregation can be found in the Appendix. 
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Gl G2 



Fig. 6 Case study 3: reaction mixtures Qi, such tliat tliey are aggregated in the 
fragment description — both states contain one protein Grb that is free on site b, one protein 
Grb that is bound to a site d of protein Sos, and one species containing a dimer of EGFR 
proteins, such that each EGFR protein is bound to one Grb protein, and one of them is 
bound to an EGF protein. Let gi S Ai C Bi and (^2 6 A2 C Bi. Then, by Theorem |9] we 
have that P{Zt = Bi) = P(Yt S {Ai,A2}) (lumpabiUty), and P{Yt = Ai) = 0.5P(Zt = Bl) 
whenever P(Yb = Ai) = P(Yo = A2) (invertability) . Moreover, by Theorem P{Yt = 
Al) — 0.5P(Zt = Bl), when t — > 00 (convergence). 



4.3 Case study 3: EGF/insulin pathway 

We take a model of the network of interplay between insulin and epidermal 
growth factor (EGF) signaling in mammalian cells from literature [5 . The 
original model suffers from the huge number of feasible multi-protein species 
and the high complexity of the related reaction networks. It contains 42956 
reactions and 2768 different molecular species, i.e. connected reaction mixtures 
which differ up to node identifiers. The reactions can be translated into a 
Kappa model of only 38 transition rules. 

The bases for the framework of site-graph-rewrite models used in this paper 
is a rule-based modeling language Kappa [7] . A Kappa rule and an example of 
the corresponding site-graph-rewrite rule are shown in Figure [5b. The general 
differences to Kappa are detailed in Section |5.5| In Figure [5p, we show the 
summary of protein interactions for this model, adapted to the site-graph- 
rewrite formalism used in this paper. Due to the independence between the 
sites a and b of protein Grb, it was proven in [5], that it is enough to track the 
copy number of 609 partially defined complexes, that are named fragments. 
Thus, the dimension of the state vector in the reduced system is 609, instead 
of 2768 in the concrete system. 

Species. Two reaction mixtures Q and Q are aggregated by relation ^iC GxG 
if they contain the same number of molecular species. 

Fragments. Let a fragment be a part of a molecular species that either does 
not contain protein Grb, or it contains only a site a of protein Grb, or it 
contains only a site b of protein Grb. Two reaction mixtures Q and Q' are 
aggregated by relation G x G if they contain the same number of 

fragments. A concrete example is demonstrated in Figure [6j 
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lumping 


# rules 


dim. 


estimated # of states 




Simple scaffold 
(3 node types) 


species 
fragment 


8 
4 


3 
2 


(n + l)(n + 2)(n + 3)/6 
(n+l)2 


Polymerization 
(2 node types) 


species 
fragment 
fragment 2 


4 
2 


n 
2 
1 


> 3P(n) 
(n+l)2 
2n+l 


EGF/insulin 
(8 node types) 


species 
fragment 


42956 
38 


2768 
609 





Table 1 Summary of the reduction for the presented case studies. In case study 1, for 
UA = TIB = iT'C = n, the number of states is reduced from 0(n'') to O(n^). The number of 



partitions of n is denoted by P{n) Ri ^ ^ 1121 . In case study 2, for ua = tlb = n, 

there is an exponential reduction in the number of states from standard to the aggregated 
CTMC. In ease study 3 (a crosstalk between the epidermal growth factor, EGF, and insulin 
pathway), the dimension of the state vector is reduced from 2768 to 609, and we did not 
estimate the size of the state space. 



5 Conclusion 

In this paper, we have studied model reduction for a Markov chain using ag- 
gregation techniques. We provided a sufficient condition for defining a CTMC 
over the aggregates, a lumpable reduction of the original one. Moreover, we 
characterized sufficient conditions for invertability, that is, when the measure 
over the original process can be recovered from that of the aggregated one. We 
also established convergence properties of the aggregated process and showed 
how lumpability and invertability depend on the initial distribution. Three 
case studies demonstrated the usefulness of the techniques discussed in the 
paper. 
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Appendix 

5.1 De-aggregation: simple scaffold 

Assume that 5 G G is such that (j)i{G) = {mAB,TnBC, rnABc)- Let niA ■= ua — 
rriAB-mABC, rriB := nA-mAB-mBC-mABC and mc := nc-mBC-niABC- 
If ^ e Aj, then au{Q) = |Aj|~^, where 



2012. 



2010. 




(10) 
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Fig. 7 Interpreting the case study 1 (simple scaffold), a) examples of reaction mixtures — 
g, Gi and 02; b) a part of the CTMC Yt £ {Ai, A2, A3, . . .}, such that g e Ai, Gi S A2, 
02 e A3; c) a part of the CTMC Zt e {Bi,B2, . . .}, such that 5 £ Bi, Si, 02 S B2. The state 
B2 is lumping of states A2 and A3. Then, by Theorem[9] we have that P{Zt = B2) = P{Yt £ 
{Ai,A2}) (lumpability), and P{Yt = Ai) = 2/3P(Zt = B2), P{Yt = A2) = l/3P(Zt = B2), 
whenever P{Yo = Ai) = 2P{Yo = A2) (invertability). Moreover, by Theorem [u] P{Yt = 
Ai) — > 2P(Yt = A2), as t — > 00 (convergence). 



The explanation is as follows. The rrij^ free nodes of type A, niB free nodes of 
type B and nic free nodes of type C can be chosen in (^'^) (J^^) (m^) possible 
ways. Among the remaining nodes, rriAB nodes of type A and itiab nodes of 
type B can be chosen in ("m"™'^) ("m"'^^) ways. There are ttiab^- different 
ways to establish bonds between ttiab identified nodes A and mAB identified 
nodes B. In the same way, we choose ttibc complexes of type (BC) among 
the hb — iriB — tuab nodes of type A, and uq — mc nodes of type B. Finally, 
there is exactly one way to choose itiabc complexes of type (ABC) among 
the UA — ruA — tuab, riB — rriB — ttiab — itibc and nc — mc — ttibc nodes of 
type A, B and C respectively. Connecting the bonds can be done in (mAscO^ 
different ways (for each node B^ , there are exactly niABC^- ways to choose the 
and ttiabc^- ways to choose C'^). The final expression follows. 

Moreover, if (f>2iG) = {mAB*,rn^Bc) and G e Mj, then a2j{G) = |Bjr\ 
where 



|B,Hf )( )m,,jf W 

\mAB*J \mAB*J ym^BcJ \m*BcJ 



We first choose the ruAB* nodes of type A and tuab* nodes of type B; 
There are m^s*! different ways to establish the bonds; In total, it makes 
(mis.) (mAB.) "^-4^*' choices. Independently, the m*BC bonds between B and 
C can be chosen in ( J^^^) („"^^)to*bc! ways. 
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( rM ^^B^ A^(p (I ) B"(a ) ( r>-( () ^ )- ^ A' ( r) 





Fig. 8 Case study 2: two-sided polymerization, a) examples of reaction mixtures; b) a 
part of the CTMC Zt S {Bi,B2, . . .}, such that Si £ Bi, 02 6 B2, c) a part of the CTMC 
e {Ci,C2, . . .}, such that Gi e Ci, 6-2 G C2. 



5.2 De-aggregation: two-sided polymerization 
Assume that s is a site-graph such that 

(l)l{s) = (Xii, . . . ,Xim,,X21, ■ ■ ■ ,X2m,X31, ■ ■ - jXsmjXil, ■ . . ,X4rn, X51, . . .,X5m)- 

We do not give the analytic expression for aii(s). For computing it, it is enough 

to use the following: 

— choosing a chain of type {A..B)i among uia nodes A and niB nodes B can 
be done in fi{mA, tub, i) = ("^^) {"^f) (*0^ ways; there are (m^ — i) nodes 
A, and (m^ — i) nodes B left. The same is used for choosing a chain of 
type {B..A)i; 

— choosing a chain of type {A..A)i among niA nodes A and tub nodes B can 

be done in ,f2{iTiA,n^B,'i) = ("i'*) ~ 1)' ways; there arc {m,A — i) 

nodes A, and {tub — {i ^ 1)) nodes B left. The same is used for choosing 
a chain of type {B..B)i; 

— choosing a chain of type {.A..B.) j among rriA nodes A and tub nodes B can 
be done in fz{mA, mB,i) = C'^'^) {'"f) (*0^/* ways; there are {niA — i) nodes 
A, and (niB — i) nodes B left. Division by i is done because of symmetries 
- every ring of type {.A..B.)i is determined by choosing i nodes of type A, 
i nodes of type B, ordering nodes A in one of i\ ways, ordering nodes B in 
one of i\ ways, but every ordering {Aji — B^i — Aj2 — -Bfe2 — ■ • • Aji — Bi~i) 
defines the same ring as {Aj2 — Bk2 — Aj^ — Bks — . . . Aji — Bki) etc. (i of 
them in total). 

Moreover, if s is such that (j)2is) = {mri,mi,a), then 



Q!2i(s) = 




If s is such that (l)2{s) = m, then 
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We choose m,.; nodes of type A among n of them, and the same num- 
ber of nodes of type B. There is niril different ways to connect them. We 
independently choose the mi,a bonds in the same way. 

To compute a3i(s), since aU of the m bonds can be either of type rriri 
or rriba, we choose i bonds of type rriri and (m — i) bonds of type mf,a, for 
i — 0, . . . ,m. 



5.3 Figure [7] 

The CTMC {Xt}, for given one node A, three nodes B and one node C contains 
different reaction mixtures over the set of nodes {A^, B^ , B^ , B^ , C^} . For ex- 
ample, let Q be the reaction mixture with the set of edges {{{A^,b), {B^, a)}}- 
There are three ways to apply the rule i?2 on Q: by embedding via function 

Vi = (^^1 Qi^: m = (y^2 (71^ ' = (^^3 c^'^^ ^ mixture with a 

set of edges {{{B"^, a), {A^,b)}, {(i?^, c), (C^, 6)}} and Q2 is a mixture with a set 
of edges {{{B\a), {A\h)}, {{B\ c), {C\h)}}, then Q{g, g^) = Q{g, g^) = c^. 
Note that = (1,0,0), 0i(gi) = (1,1,0), (j)i{g2) = (0,0,1). Let 



g e Ai, CJi e A2, £ A3. By applying the Equation (10), we have aii(tJ) = 



(tTqiMw) -1/3; "12(^1) - (iToTtTfjTw) - 1/3, and ai3 (f/2 ) - (ililolQlllO') 

1/6. 

Moreover, since 4>2{Q) = (1,0), and 02(^1) = ^2(^2) = (1, 1), let Bi,B2 G 
G2 be such that 5 sBi and ^2 £ B2. Then, a2i(g) = ((J) © l! Q QO!)-i = 

1/3 and 022(^1) = ^22(52) = ((})(i)lKl)(DlO"' = 1/9. 

Finally, observing the aggregation from Gi to G2, we have that q;i(Ai) = 

Si = 1' "2(^2) = tfti = 1/3, and a2(A3) = fffj = 2/3. 



5.4 Table [T] 

In order to illustrate how powerful the presented reduction method is in com- 
parison to the standard, species-based models, we compare the size of the state 
space in the species-based model, Gi, and in the fragment-based model, G2. 

Simple scaffold. The size of G2 is (n -I- 1)^: there are n + 1 possible situations 
between A and B nodes with 0,1,. . .,n bonds between them. The same 
holds for possible configurations between nodes of type B and C. Let f{k) 
denote the number of states with k copies of each of the nodes A, B and 
C, and with no complexes of type {ABC). If there is < i < fc complexes 
of type {AB), there can be < j < [k — i) complexes of type (BC), and we 
thus have f{k) = Yli=o(^ — i + 1) = (fc+iKfc+2) ^ rpj^^ number of complexes 
of type (ABC) can vary from to n, and thus we have the total number 
of states in Gi to be J2l=o f{k) = \ ELo(fc' + 3fc + 2) = l(ELo + 
3ELo^ + 2ELol) = 5(n(n + l)(2n + l)/6 + 3n(n-t-l)/2 + 2(n+l)) = 
(n+ l)(n-|-2)(n-|-3)/6. 
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Two-sided polymerization. We first estimate the size of G2. The value of rriri 
varies between and n, and the same holds for the value of mta- Each 
state (i, j) G {0, . . . , n} X {0, . . . , n} is reachable, since the bonds are created 
independently of each-other. The size of the state space G2 is thus (71+ 1)^. 
The size of G2 is 2n + 1, because the value of m varies between and 2n. 
Let P(n) denote the number of partitions of number n - number of ways of 
writing n as a sum of positive integers. One of the well-known asymptotics 



rii < . . . < rife, and a state Si € Gi that counts one chain of type {A..B)n^, 
one chain of type {A..B)n2 etc. It is in Gi, because it has exactly n nodes 
A and n nodes B. Therefore, the set Gi counts at least P{n) states. This 
approximation can be improved by factor three: think of the states G2 and 
Q3, which are constructed of chains of type {B..A)i, or {.A..B.)i instead of 



5.5 Relation between site-graph-rewrite rules and Kappa 

Since the main purpose of this paper is not to formally present the reduction 
procedure for a general rule-set, we described the rule-based model directly as 
a collection of site-graph-rewrite rules, which is a simplification with respect 
to standard site-graph framework of Kappa ([5]). The simplification arises in 
three aspects. 

First, the site (protein domain) in Kappa may be internal, in the sense that 
they bear an internal state encoding, for instance, post-translational modifica- 
tion of protein-residues such as phosphorylation, methylation, ubiquitylation 
- to name a few. Moreover, one site can simultaneously serve as a binding site, 
and as an internal site. We omit the possibility of having internal sites, but, 
it can be overcome: for example, the phosphorylation of a site can be encoded 
by a binding reaction to a node with a new name, for example. Ph. In order 
to mimic the standard unimolecular modification process by this bimolecular 
one, we need to ensure that the nodes of type Ph are always highly abundant, 
that is, are not rate limiting at any time. As a side remark, we point out that 
in reality it takes a binding event (e.g. binding of ATP) for a modification to 
happen. If a site is both internal and binding site, another copy of the site is 
created, so that one site bears an internal state, and another one is a binding 
state. A Kappa rule and an example of the corresponding site-graph-rcwrite 
rule are shown in Figure [5]d. 

Second, each Kappa program has a predefined signature of site types and 
agent types, where the agent type consists of a name, and a predetermined 
interface (set of sites). Each node of a 'Kappa' site-graph is assigned a unique 
name. On top of that, a type function partitions all the nodes according to 
their agent type. We instead embed the information about the node type (and 
we also abandon the use of term 'agent' in favor of 'node') directly in the name 
of the node: a node u*, i € N is of type v; The rules are accordingly written 
with these generative node names. The interface of a node type v is read 




one partition n = ni + . . . + Uf^, 



{A..B),. 
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from the collection of sitc-grapli-rewrite rules, as a union of all the sites which 
are assigned to v along the rules. Our formalism cannot specify a rule which 
operates over a connected site-graph with more than one node of a certain 
type, but the examples which we present here do not contain such rules. 

Third, we restrict to the conserved systems - only edges can be modified 
by the rules, while Kappa can specify agent birth or deletion. 

Finally, it is worth noting that wc define the notion of embedding in a 
non-standard way, through a combination of node-renaming function and sub- 
site-graph property. 



